Diversions.f90

Simulate river diversions along the river network
Source Code

!! Simulate river diversions along the river network
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.1 - 2nd September 2024  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 25/Jan/2024 | module split from reservoirs |
! | 1.1      | 02/Sep/2024 | routine DiversioSaveStatus to save discharge for  starting next simulation|
! | 1.2      | 04/Mar/2025 | no more need to assign weir-change-doy in diversion configuration file|
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! This module includes data and routines to manage 
! diversions along the river network in the distributed model.
! The module can manage both bypass channels (diverted flow 
! from a point upstream is discharged back to the same river),
! and diversion channels (diverted flow is discharged into 
! another natural drainage system nearby).
! 
MODULE Diversions

! Modules used: 

USE DataTypeSizes, ONLY : &  
  ! Imported Parameters:
  float, short

USE GeoLib, ONLY : &
  !Imported definitions:
  Coordinate, &
  !imported routines:
  DecodeEpsg

USE TableLib, ONLY : &
  !Imported definitions:
  Table, &
  !Imported routines:
  TableNew, &
  TableGetValue, &
  TableGetNrows, &
  TableSetId, &
  TableSetTitle, &
  TableSetRowCol, &
  TableSetColHeader, &
  TableSetColUnit, &
  TableFillRow, &
  TableExport

USE IniLib, ONLY: &
  !Imported derived types:
  IniList , &
  !Imported routines:
  IniOpen, &
  IniClose, &
  IniReadInt , &
  IniReadString, &
  IniReadreal, &
  KeyIsPresent

USE GridLib, ONLY : &
  !Imported definitions:
  grid_real
    
USE StringManipulation, ONLY : &
  !Imported routines:
  ToString, &
  StringTokenize, &
  StringToShort, &
  StringToFloat
  !StringCompact, StringToUpper, &
    
USE GridOperations, ONLY : &
  !Imported routines:
  GetIJ

USE DomainProperties, ONLY : &
    !imported variables:
    mask

USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

USE Chronos, ONLY : &
  !Imported definitions:
  DateTime, &
  !Imported variables:
  timeString, &
  !Imported operators:
  ASSIGNMENT (=)

IMPLICIT NONE 

! Global (i.e. public) Declarations: 

!type definition
TYPE Diversion
	INTEGER (KIND = short)  :: id !!diversion id
	CHARACTER (LEN = 100)   :: name !!diversion name
    TYPE (Coordinate)       :: xyz !!easting, northing and elevation in real world
    REAL (KIND = float)     :: xout !! x coordinate where diverted flow is discharged
    REAL (KIND = float)     :: yout !! y coordinate where diverted flow is discharged
	INTEGER (KIND = short)  :: r  !!cell row i
	INTEGER (KIND = short)  :: c  !!cell column j
    INTEGER (KIND = short)  :: rout  !!cell row  where off-stream pool outflow is discharged
	INTEGER (KIND = short)  :: cout  !!cell column where off-stream pool outflow is discharged
    INTEGER (KIND = short)  :: fileunitOut !!file unit for writing results
	REAL (KIND = float)     :: Qout !! discharge flowing in the natural river downstream diversion (m3/s)
	REAL (KIND = float)     :: eFlow (365) !! daily environmental flow [m3/s]
	TYPE (Table)            :: weir !! stream-diverted flow relationship
    INTEGER (KIND = short)  :: weirDOY (365) !!weir function used daily
    REAL (KIND = float)     :: channelLenght !! diversion channel lenght (m)
    REAL (KIND = float)     :: channelSlope !! diversion channel slope (m/m)
    REAL (KIND = float)     :: channelManning !! diversion channel Manning roughness [s m^-1/3]
    REAL (KIND = float)     :: channelWidth !! diversion channel bottom width (m)
    REAL (KIND = float)     :: channelBankSlope !! diversion channel section bank slope (deg)
    REAL (KIND = float)     :: QinChannel !! Input discharge into diversion at time t+dt (m3/s)
    REAL (KIND = float)     :: QoutChannel !! Output discharge into diversion at time t+dt (m3/s)
    REAL (KIND = float)     :: PinChannel !! Input discharge into diversion at time t (m3/s)
    REAL (KIND = float)     :: PoutChannel !! Output discharge into diversion at time t (m3/s)
    TYPE (Diversion), POINTER  :: next  !!dynamic list
END TYPE Diversion

!Public variables
INTEGER (KIND = short) :: dtDiversion
INTEGER (KIND = short) :: dtOutDiversion
TYPE (Diversion), POINTER :: diversionChannels !list of diversion channels

!Public routines
PUBLIC :: InitDiversions
PUBLIC :: OutDiversionsInit
PUBLIC :: GetDiversionById
PUBLIC :: OutDiversions
PUBLIC :: DiversionSaveStatus

!local variables
INTEGER (KIND = short), PRIVATE :: nDiversions !! total number of diversions

!Private routines
PRIVATE :: ReadWeir
PRIVATE :: SetDailyArray



!=======
CONTAINS
!=======
! Define procedures contained in this module. 

  
!==============================================================================
!| Description:
!   Initialize diversions
SUBROUTINE InitDiversions  &
!
(fileIni)

IMPLICIT NONE

!arguments with intent in:
CHARACTER (LEN=*), INTENT(IN)     :: fileIni !!diversion configuration file


!local declarations
TYPE (IniList)   :: iniDB
TYPE (Diversion), POINTER :: currentDiversion !points to current diversion
CHARACTER (LEN = 300)     :: string
INTEGER (KIND = short)    :: nArgs 
INTEGER (KIND = short)    :: k, i, j
CHARACTER (len=100), POINTER   :: args(:)
INTEGER(KIND = short), ALLOCATABLE :: doy (:)
TYPE (Table) :: initialQ ! for discharge initialization from previous simulation



!-------------------------------end of declaration-----------------------------

CALL Catch ('info', 'Diversions', 'Initializing diversion channels ')

!--------------------------------------------
!  open and read configuration file
!--------------------------------------------

CALL IniOpen (fileIni, iniDB)

!--------------------------------------------
!  allocate and populate diversion channels 
!--------------------------------------------

nDiversions =  IniReadInt ('ndiversions', iniDB)

!prepare list of reservoirs
 NULLIFY (diversionChannels)
 DO k = 1, nDiversions
   IF (.NOT. ASSOCIATED (diversionChannels) ) THEN
     ALLOCATE (diversionChannels)
     currentDiversion => diversionChannels
   ELSE
     ALLOCATE (currentDiversion % next)
     currentDiversion => currentDiversion % next
   END IF
   
   !id
   currentDiversion % id = &
       IniReadInt ('id', iniDB, section = ToString(k))
   
   !name
   currentDiversion % name = &
       IniReadString ('name', iniDB, section = ToString(k))
   
   !coordinate
   currentDiversion % xyz % easting = &
       IniReadReal ('easting', iniDB, section = ToString(k))
   currentDiversion % xyz % northing = &
       IniReadReal ('northing', iniDB, section = ToString(k))
   currentDiversion % xyz % system = &
       DecodeEpsg (IniReadInt ('epsg', iniDB))
   
   !read x and y coordinate where outflow from diversion is discharged
   currentDiversion % xout = &
       IniReadReal ('xout', iniDB, section = ToString(k))
   currentDiversion % yout = &
       IniReadReal ('yout', iniDB, section = ToString(k))
       
   !local coordinate
   CALL GetIJ ( X = currentDiversion % xyz % easting, &
                Y = currentDiversion % xyz % northing, &
                grid = mask, i = currentDiversion % r, &
                j = currentDiversion % c )
   
   CALL GetIJ ( X = currentDiversion % xout, &
                Y = currentDiversion % yout, &
                grid = mask, i = currentDiversion % rout, &
                j = currentDiversion % cout )
   
   !read weir data
   CALL ReadWeir (iniDB, k, currentDiversion )
   
   !channel lenght
   currentDiversion % channelLenght = &
           IniReadReal ('channel-lenght', iniDB, section = ToString(k) )
   !channel slope
   currentDiversion % channelSlope = &
           IniReadReal ('channel-slope', iniDB, section = ToString(k) )
   !channel roughness coefficient
   currentDiversion % channelManning = &
           IniReadReal ('channel-manning', iniDB, section = ToString(k) )
   !channel section bottom width
   currentDiversion % channelWidth = &
           IniReadReal ('section-bottom-width', iniDB, section = ToString(k) )
   !channel section bank slope
   currentDiversion % channelBankSlope = &
           IniReadReal ('section-bank-slope', iniDB, section = ToString(k) )
   
    !environmental flow
    IF ( KeyIsPresent ('e-flow', iniDB, section = ToString(k) ) ) THEN
       string = IniReadString ('e-flow', iniDB, section = ToString(k))
       currentDiversion % eFlow = SetDailyArray (string)
    ELSE !e-flow = 0.
      currentDiversion % eFlow = 0.
    END IF
          
   
 END DO

!--------------------------------------------
!  initialize discharge value 
!--------------------------------------------
 
IF ( KeyIsPresent ('path-hotstart', iniDB ) ) THEN
    !read file to initialize discharge
    string = IniReadString ('path-hotstart', iniDB )
    write(*,*) trim(string)
    CALL TableNew ( string, initialQ )
    write(*,*) 'fine TableNew'
    currentDiversion => diversionChannels
    write(*,*) nDiversions
    DO i = 1, nDiversions
        string = TRIM ( ToString ( currentDiversion % id ) )
        write(*,*) trim(string)
        CALL TableGetValue ( string, initialQ, 'id', 'Qin', &
                             currentDiversion % PinChannel  )
        CALL TableGetValue ( string, initialQ, 'id', 'Qout', &
                             currentDiversion % PoutChannel  )
        write(*,*) currentDiversion % PinChannel , currentDiversion % PoutChannel
        currentDiversion => currentDiversion % next
    END DO
END IF

!--------------------------------------------
!  close configuration file
!--------------------------------------------

CALL IniClose (iniDb)


RETURN
END SUBROUTINE InitDiversions


!==============================================================================
!| Description:
!   initialise files for output
SUBROUTINE OutDiversionsInit &
  !
  ( pathOut )

IMPLICIT NONE


!arguments with intent(in):
CHARACTER ( LEN = *), INTENT(IN) :: pathOut


!local declarations:
TYPE (Diversion), POINTER   :: current !current diversion
CHARACTER (len = 100) :: string

!------------------------------end of declarations-----------------------------

current => diversionChannels
DO WHILE (ASSOCIATED(current)) !loop trough all diversions
    
    string = ToString (current % id)
    current % fileunitOut = GetUnit ()
    OPEN(current % fileunitOut,&
    file = pathOut (1:LEN_TRIM(pathOut))//'diversion_'//&
               TRIM(string)//'.out')
    WRITE(current % fileunitOut,'(a)') 'FEST: diversion routing'
    WRITE(current % fileunitOut,'(a,a)') &
         'diversion name: ', current % name &
         (1:LEN_TRIM(current % name))
    WRITE(current % fileunitOut,'(a,i4)') 'diversion id: ', current % id
    WRITE(current % fileunitOut,*)
    WRITE(current % fileunitOut,'(a)')'data'
    WRITE(current % fileunitOut,'(a)') 'DateTime &
            Qupstream[m3/s]  Qdownstream[m3/s] &
            QinChannel[m3/s] QoutChannel[m3/s]'

     current => current % next
    
END DO

RETURN
END SUBROUTINE OutDiversionsInit
  
  
!==============================================================================
!| Description:
!   save diversion state variables on file
SUBROUTINE DiversionSaveStatus &
  !
  ( pathOut, time )

IMPLICIT NONE


!arguments with intent(in):
CHARACTER ( LEN = *), INTENT (IN) :: pathOut
TYPE (DateTime), OPTIONAL, INTENT (IN) :: time

! local declarations:
TYPE (Table) :: tab
CHARACTER (LEN = 300) :: string
CHARACTER (LEN = 100) :: row (3)
INTEGER (KIND = short) :: i
TYPE (Diversion), POINTER :: currentDiversion !points to current diversion

!------------------------------end of declarations-----------------------------

!create new table
CALL TableNew ( tab )

!populate table
string = 'diversion status'
CALL TableSetId ( tab, string)

IF ( PRESENT (time) ) THEN
  timeString = time
  string = 'diversion status at: ' // timeString
ELSE
  string = 'diversion status at the end of simulation' 
END IF
CALL TableSetTitle ( tab, string)

!Allocate variables
CALL TableSetRowCol ( tab, nDiversions, 3 ) 

!set column header and unit
CALL TableSetColHeader (tab, 1, 'id')
CALL TableSetColHeader (tab, 2, 'Qin')
CALL TableSetColHeader (tab, 3, 'Qout')

CALL TableSetColUnit (tab, 1, '-')
CALL TableSetColUnit (tab, 2, 'm3/s')
CALL TableSetColUnit (tab, 3, 'm3/s')

!fill in rows
currentDiversion => diversionChannels

DO i = 1, nDiversions
     !id
     row (1) = ToString (currentDiversion % id)
     !Qin
     row (2) = ToString (currentDiversion % QinChannel)
     !Qout
     row (3) = ToString (currentDiversion % QoutChannel)
     
     currentDiversion => currentDiversion % next
     
     CALL TableFillRow (tab, i, row)
END DO

IF (PRESENT(time)) THEN
	timeString = time
	timeString = timeString (1:19) // '_'
	timeString (14:14) = '-'
	timeString (17:17) = '-'
		
ELSE
	timeString = '                    '
END IF

string = TRIM(pathOut) // TRIM(timeString) // 'diversions.tab'

CALL TableExport (tab, string )

RETURN
END SUBROUTINE DiversionSaveStatus
  
  
  !==============================================================================
!| Description:
!   given a list of diversion channels, returns a pointer to one diversion by id
FUNCTION GetDiversionById &
( list, id ) &
RESULT (p)

IMPLICIT NONE 

!Arguments with intent in:
TYPE (Diversion), POINTER, INTENT(IN) :: list !list of reservoirs
INTEGER (KIND = short), INTENT(IN)    :: id

!Arguments with intent out:
TYPE (Diversion), POINTER :: p !pointer to reservoir

!local arguments:
LOGICAL :: found

!------------end of declaration------------------------------------------------

!loop througout list of reservoirs searching for id
p => list
found = .false.
DO WHILE (ASSOCIATED(p)) 
  IF (p % id == id) THEN
    found = .TRUE.
    EXIT
  ELSE
    p => p % next
  END IF
ENDDO

IF (.NOT. found ) THEN
  CALL Catch ('error', 'Diversions', 'diversion not found by &
                       function GetDiversionById: ', argument = ToString(id))
END IF

RETURN
END FUNCTION GetDiversionById

!==============================================================================
!| Description:
!   write results on files
SUBROUTINE OutDiversions &
  !
  (list, time, Qin, Qout)

IMPLICIT NONE


!arguments with intent(in):
TYPE (Diversion), POINTER, INTENT(IN) :: list !list of diversion channels
TYPE (DateTime),           INTENT(IN) :: time
TYPE (grid_real),          INTENT(IN) :: Qin
TYPE (grid_real),          INTENT(IN) :: Qout

!local declarations:
TYPE(Diversion), POINTER   :: current !current diversion in the list

!------------------------------end of declarations-----------------------------
timeString = time
current => list
DO WHILE (ASSOCIATED(current)) !loop trough all diversion channels

        WRITE(current % fileunitOut,'(a,4(" ",e14.7))') timeString,&
            Qin % mat(current % r, current % c),&
            Qout % mat(current % r, current % c),&
            current % QinChannel, &
            current % QoutChannel
    
     current => current % next
    
END DO

RETURN
  END SUBROUTINE OutDiversions
  
  
!==============================================================================
!| Description:
!   populate  array of daily values
FUNCTION SetDailyArray &
!
( value ) &
!
RESULT ( array )

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *) :: value


!Local declarations:
REAL (KIND = float) :: array (365)
REAL (KIND = float) :: scalar
LOGICAL :: error
TYPE (Table) :: valueTable
INTEGER :: i
REAL (KIND = float) :: doyStart, doyStop

!----------------------------end of declarations-------------------------------

!check that value is a number
scalar = StringToFloat (value, error)
IF  ( error ) THEN !value changes in time
    CALL TableNew (value, valueTable)
    array = 0.
    
    DO i = 1, TableGetNrows (valueTable)
        CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-start', &
                             match = 'exact', &
                             valueOut = doyStart )
      
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-stop', &
                             match = 'exact', &
                             valueOut = doyStop )
       
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='value', &
                             match = 'exact', &
                             valueOut = scalar )
       
       array ( INT (doyStart) : INT (doyStop) ) = scalar
        
    END DO
    
ELSE !no error, value is a scalar
    array = scalar
END IF

RETURN
END FUNCTION SetDailyArray

!==============================================================================
!| Description:
!   read weir table
SUBROUTINE ReadWeir &
  !
  (iniDB, k, div)

IMPLICIT NONE

!Arguments with intent(in):
TYPE(IniList), INTENT (IN) :: iniDB
INTEGER (KIND = short), INTENT (IN) :: k

!Arguemnts with intent (inout):
TYPE(Diversion), POINTER, INTENT(INOUT) :: div 


!local declarations
CHARACTER (LEN = 300)          :: string
INTEGER (KIND = short)         :: nDOY
INTEGER (KIND = short)         :: shortInt 
INTEGER (KIND = short)         :: i, j
INTEGER(KIND = short), ALLOCATABLE :: doy (:)
LOGICAL                       :: isString

!---------------------------end of declarations--------------------------------

string = IniReadString ('weir', iniDB, section = ToString(k))
CALL TableNew (string, div % weir)

nDOY = div % weir % noCols - 1
ALLOCATE ( doy ( nDOY ) )

j = 0
DO i = 1, div % weir % noCols
    shortInt = StringToShort ( div % weir % col ( i ) % header, isString )
    
    IF ( .NOT. isString ) THEN !number detected
        j = j + 1
        doy ( j ) = shortInt
    END IF
    
END DO

!set  when (doy day of year) geometry table changes
div % weirDOY (:) = MAXVAL (doy)
       !
DO i = 1, nDOY
    DO j = doy (i), 365
        div % weirDOY (j) = doy (i)
    END DO
END DO
       
DEALLOCATE (doy)

RETURN
END SUBROUTINE ReadWeir

END MODULE Diversions
Diversions.f90 Source File

Contents

Modules

Source Code

Source Code
